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Abstract 

We consider laminar flow in channels constrained geometrically to remain between two parallel 
planes; this geometry is typical of microchannels obtained with a single step by current microfabrication 
techniques. For pressure-driven Stokes flow in this geometry and assuming that the channel dimensions 
change slowly in the streamwise direction, we show that the velocity component perpendicular to the 
constraint plane cannot be zero unless the channel has both constant curvature and constant cross- 
sectional width. This result implies that it is, in principle, possible to design "planar mixers", i.e. 
passive mixers for channels that are constrained to lie in a flat layer using only streamwise variations of 
their in-plane dimensions. Numerical results are presented for the case of a channel with sinusoidally 
varying width. 
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I. INTRODUCTION 



The rapid development of microfluidic systems and their applications in domains such as 
aeronautics, chemistry, material synthesis, medical diagnostics and drug delivery has recently 
motivated investigations of new research questions in fluid dynamics at small scales (0,0,0]). A 
particularly active research area is the design of mixers for microdevices. This task has proven 
itself to be a real engineering challenge due to a variety of physical and practical constraints. 
Physically, the typical cross-sectional dimensions of microchannels in current microfluidic sys- 
tems are on the order of 10 — 100 /mi (0,0). On this scale, practical flows of common liquids 
(U ~ 0.1 — 10 cm/s) have low Reynolds numbers, Re < 10, and turbulent mixing of the fluid 
does not generally occur. At the same time, the Peclet number for mass transfer in these flows, 
defined as the ratio of a typical diffusive time to a typical advection time is high, Pe > 100, 
and purely diffusive mixing across the flow is slow. For applications that require mixing, it is 
therefore necessary to design channels that will lead to efficient convective mixing. 

In a channel geometry, the strongest gradients of concentration are typically oriented in 
a direction normal to the principal axis of the channel, because the gradients exist between 
co-flowing streams of distinct chemical makeup. An effective mixing flow in a channel must 
therefore stir the fluid over the cross-section with transverse flows as the fluid progresses down- 
stream in the axial flow; such a flow must possess three non-zero components. An effective 
stirring flow rapidly folds the fluid into itself so as to decrease the distance that diffusion must 
act to homogenize concentrations. The potential of a given channel-flow to mix efficiently can 
be judged by several characteristics: (1) the ratio of the transverse to the axial velocities. If this 
ratio is too small, then the axial length of channel required for mixing is likely to be imprac- 
tically long, regardless of the detailed character of the flow; (2) the distribution of transverse 
flows within the cross-section of the channel. If the flows are confined to small areas within 
the cross-section, then they will be ineffective at exchanging fluid between these areas and with 
areas in which no flow exists. (3) the evolution of the transverse flows as a function of axial 
position along the channel. An efficient mixer should produce Lagrangian chaos; in particular, 
no streamwise symmetries should be present . This feature is achieved in a channel-geometry 
when the position and orientation of transverse flows vary axially such that all fluid elements 
travelling in the channel experience an alternating sequence of rotational and extensional flows 
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In designing a laminar mixer for microfluidic applications, it is also important to take into 
account the constraints imposed by technology and conversation laws. Designs of mixers should 
be scalable to smaller systems, resistant to fouling by particulate matter, and efficient with 
respect to power consumption. Moreover, the typical lithographic processes currently used in 
microfabrication lead to Dlanar geometries. 

A variety of active (Q, 0) and passive methods have recently been 



proposed to generate stirring flows for mixing purposes. All of the passive designs that have 
been demonstrated to be effective rely however on multilayer or non-planar structures in order 
to generate three-dimensional flows. A single lithographic step generates a single flat layer of 
structure with horizontal top and bottom walls and (relatively) vertical side walls (see Figure 
([14|). More complicated structures, e.g. a non-intersecting cross-over between two channels, 
require multiple lithographic steps, with spatial alignment between each layer of structure. Each 
added layer of structure increases the difficulty of fabrication and complicates scaling down to 
smaller devices. 

Motivated by such practical considerations in this paper, we will be interested in character- 
izing the potential for mixing of the simplest planar geometrical configurations obtained in a 
single fabrication step. Two such geometries are illustrated in Figure ^ Suppose we fabricate 
a channel constrained between two parallel planes of constant separation, can such a configu- 
ration mix? In order to be able to give an answer to this question that is robust to change in 
flow conditions, we will assume zero Reynolds number flows in the channel. This assumption 
also implies that our conclusions will remain valid in smaller flow systems of the same design. 

It is known that the steady-state velocity field in a straight channel of constant rectangular 
cross section is unidirectional (^^) and therefore cannot mix except by molecular diffusion; 
similarly, the velocity field in a curved channel of constant cross-section and constant curvature 
is unidirectional ( lg])0|- As a consequence, the simplest potential design for a steady-state 
mixer for Stokes flow is that of a channel with variations of shape, that include changes of both 
curvature and cross-sectional dimensions in the streamwise direction. 

In order for flow in such a channel to potentially mix by advection in the three dimensions 
of the channel, the velocity field needs three non-zero components. While it is clear that these 
variations in shape will lead to non-zero in-plane components of the velocity, as would also 
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be the case in a truly two-dimensional channel, it is not obvious that the (third) out-of-plane 
component will always be non-zero. We ask therefore the following question: under which 
circumstances is the out-of-plane component of the velocity field always non-zero? And in this 
case, what is the expected magnitude of the vertical flow? 

The flows in a circular pipe of varying cross-section (^J) or varying small curvature 
have been studied and three-dimensional flow is obtained at zero Reynolds number. However, 
because the equation for the shape of a circular pipe couples the two directions that are perpen- 
dicular to its axis of symmetry, these results cannot be applied to the flow in a planar geometry 
and a separate analysis has to be carried out. Recently, Balsa studied the secondary flow 
in a Hele-Shaw cell in which a vertical cylinder is immersed, at Reynolds number unity based 
on the cylinder length, and showed the presence of streamwise vorticity in a boundary layer on 
the cylinder surface; an earlier study by Thompson j^J focused on viscous features. 

The geometry of a generic microchannel constrained between two parallel planes with fixed 
separation and with no obstacles is illustrated in Figure ^ (top). The shape of the channel can 
be entirely described by two degrees of freedom: (1) the shape of its centerline plane and (2) 
the local symmetric width of the channel around this centerline. We will consider in this paper 
the consequences of both and will treat each of them separately for simplicity. 

The paper is organized as follows. In section Hi Al we consider the case of a straight channel 
with varying cross section in the direction perpendicular to both the flow and the constraint 
plane and in section III Bl we consider the case of a curved channel of constant cross section 
but varying curvature. In both cases, under the assumption of an arbitrary but slowly varying 
cross section and curvature respectively, we show that the velocity component perpendicular 
to the constraint plane cannot be zero unless cross-section and curvature are both constant, 
and therefore the flow is fully three-dimensional in all other cases. We apply these results in 
section HTT1 where we calculate the leading-order velocity field in the case of a straight channel 
of varying cross section and illustrate the flow patterns on a sinusoidally varying channel. We 
conclude in section IIVI with a discussion of both the practical advantages and limitations that 
these results imply for mixing design. Appendices IA II and IA2I present proofs for some of the 
results used in sections fll Al and III Bl respectively and Appendix IB1 presents the calculation for 
the leading-order velocity field in the case of a channel of arbitrary shape. 
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II. THREE-DIMENSIONALITY OF THE FLOW 



A. Straight microchannel of varying cross section 

In this section we consider the case of a straight microchannel of varying cross section, as 
illustrated in Figure El A pressure-driven flow takes place in the x direction of a channel of 
constant height 2h in the ^-direction and varying width 2hf(x/\) in the y direction, where A 
is the axial length scale on which such variations occur. The equations of motion and mass 
conservation for an incompressible Stokes flow are written 

Vp = /iV 2 u, V.u = (1) 

with the no-slip boundary condition u(x, y = ±hf(x/X), z) = u(x, y,z = ±h) = 0, on the four 
bounding surfaces. The flow rate Q is set by upstream conditions and is constant 

rh p+hf(x/X) 



u.e x )dydz = Q. (2) 

h J-hf(x/X) 

We now make the assumption that the width of the channel is slowly varying. If we define 
e = h/X, this assumption is equivalent to assuming that e C 1. Using the notations u = (u, v, w) 
for the velocity field, we can now nondimensionalize the previous set of equations CQ)-(J2I) by 
scaling lengths, velocities and pressure by 

(?,y,z) = (Xx,hy,hz), (u,v, w) = j^(u, ev, ew), p = ^rP- (3) 

Dropping the tildes in the dimensionless variables for convenience, the dimensionless Stokes 
equation is 

/ „ r) 2 f) 2 & \ r _ - i 

(4) 



u, e 2 v , e 2 w 



' \ dx 2 dy 2 dz 2 

and the dimensionless continuity equation is 

du dv dw . . 

dx dy dz 

We look for a regular perturbation expansion for both the velocity and pressure fields under 
the form 

(u,v,w,p) = (u ,v ,w ,po) + e 2 (u 2 ,V2,w 2 ,p2) + C(e 4 ) (6) 
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which is usual in lubrication theory (see e.g. [2JJ] ) . The leading-order O(e ) term of Stokes 
equation (JTJ) is given by 

dp ( d 2 d 2 \ dp dp 



dx \dy 2 dz 2 J dy dz 

subject to the no-slip boundary condition uo(x,y = ±f(x),z) = uo(x,y,z = ±1) = 0, and 
constant flow rate 

r 1 rf( x ) 

4/ / u dydz=l. (8) 



JO JO 

The axial velocity u is then easily calculated by separation of variables (|l5j) 

u (x,y,z) = -—^z -l) + ^^ r [ coshknf{x) )cosk nZ ^ (9) 

with k n = (n + |)7r. This solution is simply the axial Poiseuille velocity in a straight channel 
of constant cross section evaluated at each location along the channel. The leading-order axial 
pressure gradient is then given by the flow rate condition (jSJ) which leads to 

dpo 3 I 6 taxih(k n f (x)) 



dx Af(x){f(x)^ o k\ 

Note that (jJJ) and © show that at next order in e 2 , the leading-order out of plane velocity 
field (fo,wo) satisfies a two-dimensional Stokes equation with an effective distribution of mass 
sources and sinks given by — dito/dx. Let us now assume these sources and sinks lead to a 
planar velocity field is planar in the sense that the component of the velocity perpendicular to 
the constraint plane is zero, wq = 0. In this case, the continuity equation from (0) is written 

^ + ^ = 0, (11) 
ox ay 

and allows us to solve explicitly for the ^/-component of the velocity Vq. Using the fact that u 
satisfies the no-slip condition on the walls of the channels, it is straightforward to obtain that 
v o is given by 

vo{x ) y ) z) = —l u (x,y',z)dy'. (12) 

dx J-m 

The solution (j!2)l satisfies the no-slip conditions for vq at z — ±1 and y = —f(x). If it also 
satisfies the remaining no-slip condition at y = f(x) then the leading-order solution would be 
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entirely characterized and the flow would be planar at leading-order. The condition at y = f(x) 
will however be satisfied if and only if 

d f f{x) 

— / u (x,y ,z)dy = (13) 

OX J-f(x) 

for all values of x and z in the channel. Using the solution J5J), (|13|) can be integrated once to 
get 

■(tanh(knf(x)) - k n ) cosk n z } = (14) 




A;: 1 



n 



In order for equation (pPfj) to be satisfied for all \z\ < 1 and x > 0, it is then necessary that for 
all n > 0, 

^-(tanh(k n f(x)) - k n ) = 5 n (15) 
da; 

where the {6 n } are constants independent of x. As is shown in Appendix IA 11 this can only 
be true if f(x) is constant, i.e. if the channel cross section is constant. As a consequence, 
the vertical component of the velocity field u>o cannot be zero unless the cross section of the 
channel is constant, in which case vq = wq = 0. When this is not the case and the cross section 
is changing along the channel, then the two-dimensional solution (J§|)- (fT2|) is inconsistent and 
the velocity field is fully three-dimensional at this order. 



B. Channel of constant cross section with varying curvature 

We now proceed in the same manner as in section III Al for the case of the channel of constant 
cross section but varying curvature, as illustrated in Figure El A pressure-driven flow takes 
place in the axial direction, denoted as s, of a channel of constant height 2h in the z-direction 
and constant width 2d in the third direction, denoted as n for "normal". The centerline of 
the channel is not straight but curved with local radius of curvature R{s) = Rof(s/\) in the 
orthogonal (e n , e s , e z ) frame, where A is the typical length scale along the channel on which this 
local curvature changes. In this geometry, and using the notations u = (u, v, w) for the velocity 
field, it is possible after some algebra to write the dimensional Stokes equation under the 
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form 



1 R(s) dp d 2 u ( R(s) \ 2 d 2 u d ( 1 d /rt , , 



/i R(s) + n ds dz 2 \R(s) + n J ds 2 dn \R(s) + ndn 

2R(s) dv R(s) dvdR 
+ (R{s) + n) 2 ds + {R{s)+ny^^^ , ( ' 

I dp d 2 v ( R(a) \ 2 d 2 v d ( 1 d , n , s 

+ [ n , < ^ + r LM — (R(s) + n)v 



fidn dz 2 \R(s) + n J ds 2 dn \R(s)+ndn 

2R(s) du R(s) dR ( dv . 

>n— -u , (16b) 



(R(s) + n) 2 ds (R(s) + nf ds \ ds 
1 dp d 2 w ( R(s) \ 2 d 2 w 1 d s . <9wA 



/zcfe <9z 2 \i?(s)+ny (9s 2 R(s)+ndn \ dn J 

nR(s) dRdw 

+ {R{s) + nf^s"ds~ , ( 



and the continuity equation is written 



no V % i_a ((fl( , ) + B) „ ) + ft£ = . (17) 



R(s) + n J ds R(s)+ndn dz 

The two sets of equations ()16|) and ()17|l are associated with the no-slip boundary condition on 
the walls of the channel u(s, n = ±d, z) = u(s, n, z — ±h) = 0, as well as with the condition of 
constant flow rate along the channel 

-h r d 



{u.e s )dzdn = Q. (18) 

-h J-d 

As in section Hi A| we now assume a slowly varying curvature, i.e. we assume that both ?i/A < 1 
and d/X -C 1. Equations (fTHj) - (fT£J) can be nondimensionalized by scaling lengths, velocities 
and pressure by 

(s,n,z) ~ (\3,dn,hz), (u, v, w) = ^-(u, -v, ew), p = j^jP (19) 

where we denoted the aspect ratio a = h/d = 0(1). Defining (3 = d/Ro and e = h/X <C 1, 
and dropping the tildes in the dimensionless variables for convenience, the dimensionless Stokes 



8 



equation is 



f(s) dp d 2 u i 2 ( f(s) \ d 2 u d ( a 2 d 



f(s)+(3nds dz 1 ^^ \f(s)+(3n) ds 2 + dn \f(s) + (3n dn^^ + ^ U 
2e 2 Pf(s) dv e 2 f(s)f'(s) dv 



(/(a) + /3n) 2 9s (/(s) + /3n) 2 9n ' 
dp £^ 6 4 / f(s) \ 2 8 2 v 2 d_f 1 g_ \ 

* + «f„i* „), (20b) 



(/(s) + f3n) 2 ds a(f(s) + (3n) 3 \ ads 
dp 2 d 2 w 4 / f(s) \ d 2 w e 2 a 2 d ( dw 
dz ~ 6 lh 2+e \f(s)+Pn) ~ds T + f(s) + (3ndn~ \}^ 8) + Pn) dn 
ne^f(s)f(s) dw 

(f(s)+pn) 3 ds' 

and the dimensionless continuity equation 



(20c) 



f(s) du 1 d , „, , „ , dw , , 

+ T7 - — 77- + — = 0. (21) 



f(s) + /3nds^ f(s) + (3n dn U 1 ; "~ ^ > "~ dz 
Note that j3 is not necessary small in actual MEMS applications 22j . We then look for a regular 
perturbation expansion for both the dimensionless velocity and pressure fields under the form 

(u,v,w,p) = (u ,v ,w ,p ) + e 2 (u2,V2,W2,p 2 ) + C(e 4 )- (22) 



The leading-order O(e ) of the Stokes equation (JTfiJl is 

f(s) dp d 2 u d ( a 2 d 



(f(s)+(3n)u ) (23a) 



f(s) + (3n ds dz 2 dn \ f(s) + (3ndn 

< 23b > 

together with the no-slip boundary condition uo(s,n = ±l,z) = uo(s,n,z = ±1) = 0, and 
constant flow rate 

/ u dzd7i = l. (24) 
-l J-i 

Using separation of variables, it is possible to solve for the axial velocity component in equations 
(l2*3J) - (j2H) , similarly to what was done by Rieger & Sestak (1973) for the case of a curved 
rectangular channel of constant curvature. We obtain 



dp J z 2 -l ^ fk n (f(s)+Pn) . ■ ; 

«o(s, n, z) = f(s) — <j 2{m+(3n) + |> { ) cos knz } , (25) 
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where the set of functions U n are defined by 

U n ( V ) = E n {s)K x ( V ) + F n (s)h ( V ) , (26) 

with Ii and K\ as the order-one modified Bessel functions of the first and second kind respec- 
tively, and with 

En(s) = {(f(s)-^I 1 (k-(s))-(f(s) + ^)I 1 (kt(s))}G n (s), (27a) 

Fn(s) = {(f(s) + P)K 1 (k+(s))-(f(s)-^)K 1 (k-(s))}G n (s) (27b) 

where 

= fc3(/( 8 )2 _ j g2 ) {hiKis^k+is)) - hik+is^K^k-is))}- 1 , (28) 

and 

± = (29) 

Using the identities K' = —K\ and I' Q = —I\ and the flow rate condition (|24|) we obtain the 
pressure gradient 

with 

H n {s) = E n (s) {K (k-(s)) - K (k+(s))} + F n {I (k+(s)) - I (K(s))} • (31) 



As in section III Al let us now make the assumption that the flow is planar, i.e. that the 
leading-order vertical component of the velocity field is zero, wq = 0. In this case the continuity 
equation (JT7|) becomes 

d -£ + W)l {f(s)+pn)v °=°- (32) 

which can be used to solved exactly for t> 

v (s, n,z) — — ^- < I u (t,n' , z) dn'X ■ (33) 



f(s)+i3nds u _! j 
The solution (J33)) satisfies the no-slip boundary condition at z — ±1 and n = — 1; if the 
condition at n = 1 was also satisfied, the leading-order velocity field would be two-dimensional 
at leading-order, wq = 0. The no-slip boundary condition evaluated at n = 1 will however be 
satisfied if and only if 

^u o {t,n',z)dn'\ = 0, (34) 
10 



for all values of s and z. Using the solution for the axial velocity (J25j) . (|34|) can be integrated 
once to obtain 

nM^kJ^ — ^ll_)L* w , (35 ) 

where H n is defined in (|3~TJ) . In order for (f33j) to be satisfied for all \z\ < 1 and s > 0, it is 
necessary that all for all n > 0, 

d Po { aH n (s) 2(-l)» . 

where the {7 n } are constants independent of s. As is shown in Appendix IA 21 (|3*Bj) can be 
satisfied if and only if f(s) is constant, i.e. if the curvature of the channel is constant, in 
which case v = wq = 0. When this is not the case and the curvature is changing along the 
channel, then the two-dimensional solution (|25 j) - (|33j) is inconsistent and the velocity field is 
fully three-dimensional at this order. 



III. ILLUSTRATION OF THE THREE-DIMENSIONAL FLOWS 

We have demonstrated in the previous two sections that flows in channels constrained to 
remain in a layer of constant thickness are in general three-dimensional, i.e. they possess a 
non-zero component of the velocity perpendicular to the constraint plane. We illustrate these 
results in this section for the case studied in section Hi Al of a straight channel of varying width. 
We calculate the three components of the leading-order velocity field (u ,v ,Wo) and illustrate 
the flow patterns in a sinusoidally varying channel; the calculation for the general case of an 
asymmetric channel is more intricate and we present it in Appendix [Blfor the interested reader. 



A. Governing equations 

Because the velocity field (uo,vo,wq) is three-dimensional, the continuity equation in (Jljl 
becomes 

du Q Ovn dlVn . . 

7T + 7T + ~fT = °> 37 
ox ay oz 
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where Uq is still given by equation (J0J). Under Stokes flow conditions, the pressure is harmonic 
V 2 p = 0, and therefore the velocity field always satisfies the biharmonic equation V 4 u = 0. 
Consequently within the lubrication approximation (j3J) and (JHJ), the y and z component of the 
dimensionless velocity field satisfy 

W 2 + a* J (38a) 

Vl^o = 0. (38b) 

Similarly, the vorticity is harmonic V 2 uj = 0, so that under the lubrication approximation, its 
leading-order axial component u = — ^ satisfies 

V> = 0. (39) 

It is necessary to solve the set of equations (p?7|) , (jHEJ) and (jHEJ) along with the no-slip boundary 
conditions in order to obtain the final solution for (t> , Wq)- 

B. Subset of equations 

Let us now show that it is sufficient to solve equations (}3*Tj) and ()38a|) to obtain (j38b|) and 
(J3n|l . Let us suppose (fHTj) and ()38aj) are satisfied. Evaluating the bi-harmonic of the 
continuity equation (|57j) leads to 

^Viw = -> Vl^o = T(x, y). (40) 

Because of the symmetries in the configuration illustrated in Figure ^ the solution of Stokes 
equation has to satisfy w (x,y, —z) = —w (x,y, z) (and also v (x,y, —z) = v (x,y, z)). Conse- 
quently V\wq is also odd with respect to z and necessarily F(x, y) = 0, so that equation (J38bj) 
is satisfied. In the same fashion, it is straightforward from (j3*Tj) to obtain 

-^ = 0^uj = A(x,y). (41) 

Using the fact that wo and Vq are respectively odd and even with respect to z, it is clear that u; 
is odd with respect to z so that A(x, y) = 0. The result of equation (|3T?J) is therefore recovered. 
As a consequence, it is sufficient to solve equations (|37j) and (j38aj) to obtain the complete 
solution for the leading-order velocity field. 
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C. Velocity field calculation 



To obtain the leading-order solution for the three-dimensional velocity field in the channel, 
we solve equations (|38ajl and (jSTj), along with the no-slip boundary conditions for t> and wo 
and with the axial velocity Uq given by Q. In order to do so, we use the technique introduced 
more than a century ago by Lame (|23j) to solve planar elasticity problem where biharmonic 
equations arise (see also the general discussion in Here we effectively demonstrate that 

these ideas also apply as well to slowly varying flows. Using the following symmetries in the 
velocity field 

Vo(x,y, -z) = v (x,y,z), v (x, -y, z) = -v (x, y, z), (42a) 
w (x, -y,z) = w Q (x,y,z), w Q (x,y,-z) = -w (x,y,z), (42b) 



we look for a solution of ()38a|) under the form of a double Fourier series in y and z 

v (x,y, z) = y^A n (a;, y) cosk n z + ^ B m (x, z) sin (jt^J , ( 43 ) 

n>0 m>0 \J\ X )/ 

with l m = mir and k n = (n + 1/2)tt. In order for (J38a|) and (|42aj) to be satisfied, the functions 
A n and B m are given by [28] 

A n (x, y) = a n (x)P n (x, y), B m (x, z) = b m (x)Q m (x, z) (44) 

with 

P n (x, y) = f(x) smh(k n y) - y cosh(k n y) twah(k n f(x)) (45a) 

Q m (x,z) = tanh ( 4t\] cosh ( TTT J ~ 2simi ( TTT J ( 45b ) 

and where both (a n (x)) and (b m (x)) are unknown functions to be determined. With the axial 
solution Q written 

u (x,y,z)="£u n (x,y)cosk n z, U n (x,y) = ( co^/(x)) " *) (46) 

and with (jlHJ), integration of the continuity equation (J37|) leads to the third component of the 
velocity field 

1 fdU n , , ,dP n \ . , 

\ ^ @-mbm{.%)T m (x } z) ( ( m y \ , . 
w (x, y , *) = - ^ _ ^_ + a n (x)— j sm k n z - ^ ^ cos j (47) 
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with 

- (f> ^ + If) ... (^) - ffiU* ( ) • («, 

The sets of unknown functions (a n ) and (6 m ) are finally determined by enforcing the no-slip 
boundary condition for the solution (}4Tj) at both z = ±1 and y = ±f(x). As is usually the 
case for problems involving biharmonic equations (e.g. [2^] ) , the final result involves an infinite 
system of linear algebraic equations given by 

) + 2J Aim(^m(s) (n > 0), 
m>0 

fcm(^) = fem(^) + S ^B mn (x)a n (x) (m > 0), 



(49) 



with 



n>0 



„ (dp n 



On(x) = — — (x,f(x)) i—(x,f(x)) j , (50a) 
4wn(a0 = 2< ~ 1 ^ a .^ mfcw j"^^'-^)} / Tm ( x ' 2; ) sin^n^d^;, (50b) 

Ux) = im^) g -is— /„ ^ y) cos ItmJ dy ' (50c) 

2(-l)"+ 1 /-/flP B f4»v\, 
B mn {x) = — — / -= — [x,y cos — - dy. 50d 
k n £ m T m (x, 1) J dy \f(x)J 

Note that, at a given position x, both the (a n (x)) and (6 n (a;)) are entirely determined by the 
instantaneous values of f(x) and f'(x); each subsequent order in the long- wavelength expansion 
(JHJ) will bring an additional dependance on a higher derivative of f(x). 

D. Further calculations 

1. Axial vorticity 

Given the sets of {a n } and {b n }, we can evaluate the axial component of the vorticity 

_ dwp _ dvp . 



n>0 



g Y (jix&y + an ^~dy^ ~ ^n^Pnix.y)^ sin/c„z (51a) 
■ \ ^ ( £mbm(x)T m (x, z) , ,dQ m \ . f £ m y \ 

+ £( — jw bm{x) -dr) 8m {wv' (51b) 



m>0 
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2. Quadrant-averaged velocities 



The quadrant-averaged velocities can also be evaluated, for example in the quadrant (y > 
0, z > 0). The flow rate condition (jHJ) leads to a constant average axial velocity, < uq > (x) = 
1/4. Integration of and (}4Tj) across the quadrant leads to 

^ (-!)"«„. /-/W &m/ ( x )(i + (_i^ +1 ) yi 
< v > (x) = / P n {x,y)dy + > « / <2m(x, *)dz, (52) 



n>0 J v y " u m>0 

and 



< w > (x) = 



2(_i)»+i d 
u /,l/'(.r) d.r 



^ (tanh(A; n /(x)) - k n f(x)) 



(53) 



E. Case of a sinusoidally varying wall 

We chose to illustrate the flow patterns in the case where the wall shape is described by 
the dimensionless function f(x) = 1 + 0.7 sin x; recall that the actual dimensional wall shape 
is described by /(ex) where e = h/X. The infinite system of linear equations (j4H|l was solved 
numerically by truncating it at finite values of n and m. The integrals in ()50|) only involve 
linear and trigonometric functions and are evaluated exactly. Note that apart from the system 
f!49|) . summations are also involved in equations (fTU|) for the pressure gradient and ()50c|) for 
b m (x), and they also require numerical truncations. 

For each case, numerical results were obtained, the truncation was refined and the results 
were found to converge quickly to a final solution. The truncations at n = 50 in equations 
()10|) and ()50c|) were found to be suitable to obtain the final solution. Further, a truncation at 
n = m = 20 in the infinite set of linear equations equation (}4T?|) was also found to be appropriate 
to resolve the flow fields, with results essentially unchanged for higher truncation numbers. 

Such techniques allow us to obtain everywhere in space the three leading-order velocity 
components and therefore, with a simple time advancement scheme, to follow the motion of 
individual fluid elements and obtain streamlines. 

The main results of our flow illustration are displayed in Figures HI El and El FigurelUpresents 
in-plane velocity plots at three locations along the channel direction, as well as iso- value maps 
at these locations for both in-plane velocity components (vo,Wo) and for the axial component 
of the vorticity (o>o) [2^. Figure El displays the flow streamlines along the expansion part of the 
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channel (3ir/2 < x < 5tt/2). Finally, Figure El displays the maximum cross-sectional as well as 
average value of the three components (uq,vq,wo) of the leading-order dimensionless velocity 
as a function of the location along the channel centerline. 

The numerical results confirm that the flow at leading-order is fully three-dimensional. The 
plots in Figure 0] allow us to visualize the regions of high and low velocity and vorticity and 
the streamlines in Figure El show the fluid elements are indeed vertically displaced as they are 
advected along the channel. Note that the similar plots for the contracting part of the channel 
were not included here as they can be deduced from those in Figure 0] and E] by symmetry of 
Stokes's equation. 

We also note in FigureEJthat the qualitative picture for the iso- values of v do not vary much 
between the point of minimum width (x = 3tt/2) and the point of maximum width (x = 5ir/2). 
In contrast to Vq, the distribution of vertical velocity wq is modified appreciably: it changes 
from a monotonic variation across the channel (left picture in Figure Efc) to a variation with 
local minimum/maximum in the middle of the channel and global maximum/minimum near 
the channel walls (middle and right picture in Figure Eh). Moreover, as can be seen in Figure 
the axial vorticity is maximum at the top and bottom walls and decays towards the middle 
of the channel (z=0); the contracting part is therefore the position along the channel where the 
strongest stirring of material surfaces would occur. 

Further, the results of Figure |B1 show that under the lubrication approximation, the 
magnitude of the vertical flow component Wq decreases monotonically during an expansion 
(37r/2 < x < 57r/2); by symmetry of Stokes's equation, wq increases in a similar fashion during 
a contraction of the channel, (tt/2 < x < 37r/2). 

We see also that for the particular case considered here, and within the lubrication ap- 
proximation, the leading-order y-component of the velocity vq is always about one order of 
magnitude smaller that the axial component Uq and that the vertical component Wq is about 
one order of magnitude smaller than t> ; back to the dimensional variables, these statements 
become v ~ eu/10 and w ~ eu/100. 

Finally, the integrated effect of the vertical flow along the channel length is illustrated in 
FigureElby the vertical deflection of streamlines. The deflection is larger far from the horizontal 
centerplane (see Figure Eb) and far from the vertical centerplane (see Figure Et)- Over the 
channel half period, a fluid element on the upper right streamline in Figure Ek experiences a 
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vertical displacement of about 10% of the channel half-height. 



IV. CONCLUSION 

We have shown in this paper that the only planar channel shapes for which the velocity field 
is two-dimensional under Stokes flow conditions have both constant curvature and constant 
cross section (in which case the flow field is in fact unidirectional). In all other cases for 
the variation of the cross-section and curvature, the velocity is fully three-dimensional at zero 
Reynolds number and could in principle be used to mix species in simple microdevices that can 
be manufactured with one step of microfabrication. 

A qualitative summary for the occurrence of the third component of the flow can be given 
using the two- dimensionality condition, i.e. equation ()13|) or ([34)1 . The velocity field remains 
two-dimensional in the channel if the two-dimensional flow rate Q(z) = J udy is constant along 
the channel for each z. When this is not the case and Q(z) is streamwise-dependent, a vertical 
velocity component is induced by mass conservation. What our study shows is that, under 
the lubrication approximation, the only channel geometries with no embedded obstacles for 
which this is not the case are those of both constant width and constant curvature. Note that 
alternatively, the presence of obstacles such as cylinders in an otherwise straight channel would 
provide similar geometric features necessary for the occurrence of a three-dimensional flow jsoj ] . 

As the general form of the continuity equation shows, the magnitude of the ratio of the out- 
of-plane velocity component w to the axial component u scales as the ratio of the cross-sectional 
length scale h to the length scale A over which the variations of the channel geometry occur, 

W k frA\ 

— « - = e. (54) 
u A 

The numerical results presented in section IIIIEI for a sinusoidal change in cross-section show 
that the prefactors for this scaling is about 0.01 for the ratio of the leading- order velocity fields 
Wo to Uq and indicate poor mixing. For the case h ~ A, we could expect however all orders in 
the perturbation expansion to contribute in a non trivial way, and we expect therefore that with 
this simple design a vertical flow of strength comparable to the axial flow could exist; if that 
is not the case and the prefactors for the full calculation are not of order one, the channel will 
likely present poor mixing characteristics. Note that as the Reynolds number in micromixers 
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is not exactly zero but can be as high as 100, we also expect in this case the occurrence of 
non-trivial Dean flow-like contributions to the vertical flow. 

We propose to design " planar mixers" by a succession of n mixing cells of length A along a 
single channel. In each cell, we expect the integrated displacements 5y and 5z of fluid elements 
in the cross-section, advected by the flow at velocity ?7 ax i a i, to be given by 

5y « 5z « t±U±, (55) 

where t± is the residence time for the flow in the cell t± ~ A/L^ai and U± is the magnitude 
of the transverse flow, at most U± ~ hU^^/X so that 8y ~ Sz ~ h. Since the total length 
of the mixer is nX and the displacements 5y and 5z are independent of the cell length, small 
cells A ~ h should be chosen. The challenge in the mixing design would then concern (1) the 
design of each cell, i.e. the variations of its radius of curvature and its cross-section, in order to 
obtain the maximum cross-sectional displacement and (2) the setup of the cell succession in a 
way that mixing adds up instead of cancelling out; for example the channel studied in section 
IIII El would obviously make very poor mixing cells because by symmetry of Stokes flow, every 
fluid stirring taking place in one part of the channel would be unstirred in the other part of the 
channel located immediately downstream. In general, good performance may be achieved by 
avoiding any geometrical symmetry along the streamwise direction. 

The calculations presented in this paper assumed slowly varying cross-sectional and curva- 
ture change along the channel, e <C 1. As was shown by Lucas (1972) for two-dimensional 
channels of varying shape, the regular perturbation expansions (jOJ or (|2*2*|) are expected to have 
have order one or larger radius of convergence in e; as a consequence, the conclusions reached 
using the leading-order velocity fields are valid for the entire velocity field when e is 0(1), and 
presumably higher even though our results cannot be applied directly. With current microfab- 
rication techniques, the minimum in-plane dimension (« A) that can be generated is typically 
greater than the minimum out-of-plane dimension (« h). As a consequence, the cross-sectional 
dimensions of microchannels tend to satisfy the criterion, e < 1, and the results obtained in 
this paper are expected to apply for all such cases. 

The limitation of the passive mixing strategy proposed here lies in the top-bottom symmetry 
for the velocity field, (u, v, w)(x, y, — z) = (u,v, —w)(x,y, z), due to the symmetries of the 
Stokes equations. Mixing can therefore not be achieved between the fluids located in the z > 
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and z < planes and consequently, the streams of solutions that are to be mixed must be 
introduced at the inlet of the channel with alignment in the vertical direction. The case of a 
straight channel of varying section studied in section Hi Al also possess a right/left symmetry, 
(u,v,w)(x, —y, z) = (u, —v,w)(x,y, z), (see Figure EJ) and therefore cannot mix species fluids 
located in the y > and y < planes. The configuration studied in section III Bl does not 
possess such a symmetry and should be used to transport fluid between the n > and n < 
planes, similarly to what was achieved in Stroock et al. (2002). A fully numerical approach 
to the problem (using e.g. a boundary element method or a commercial code) would allow a 
detailed study of the proposed mixing design, its optimization, and dispersion characteristics. 
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APPENDIX A: PROOFS 

1. Proof of result (fT3)l for straight channels 

We show in this section that the only set of width functions f(x) that satisfy (|15|) are the 
constant functions. Let us assume that is satisfied for a function f(x). We first rewrite 
equation (fTUj) for the pressure gradient under the form 



Substituting the expression obtained in (|A2|) into (jAljl leads to a closed form solution for the 
axial pressure gradient 




(Al) 




(A2) 




(A3) 
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As a conclusion, the functional form (jA3|) obtained for the pressure gradient is not consistent 
with that given by the assumption of the two-dimensionality of the flow (fT5)l 

dp 5 n 



(A4) 



da; tanh(/c n /(x)) — fc n ' 
unless the function f(x) is constant. 

2. Proof of result for curved channels 

We show in this section that the only set of curvature functions f(s) that satisfy (|36j) are 
the constant functions. Let us assume that (jHfijl is satisfied for a function f(s). We first rewrite 
equation (JHUj) for the pressure gradient under the form 

^t^-^-WMw^l)-- 1 - (A5) 



We then rewrite the condition ()3bj) for two- dimensionality of the flow as 

dp aH n (s) 2(-l) n dp 
f s 1 1 = 77, g "7T~ + A6 



,/(«)-/3, 

Substituting (|A6|) into (|A5|I leads to a closed-form solution for the streamwise pressure gradient 
f(s)^= A3 !^ 7 T" / , A3 = 3/3<J>;^^ Ii -^ ) A 4 = 6>;— (A7) 



ds In 2 



77w+/A a - ' A3 " 3/5 \^^—~ "2 ' A4 ~ 6 ^' 

^ J ~ A 4 U>0 n J n>0 n 



further, it is possible to use the asymptotic behaviors near x ~ oo of Bessel functions (see e.g. 

26]) I p (x) ~ e x / \ / 2ttx, K q (x) ~ ire~ x /^/2x to obtain the asymptotic behaviors of E n (s), F n (s), 
G n (s) as n — > +oo, from equations (|27a|l . (J27bj) and (f28j) respectively. It is then straightfoward 
to obtain the asymptotic behavior of H n 

4(-l)"/[» 

H " (s) (A8) 

This behavior and the condition of two-dimensionality (J36j) allows to obtain an alternate func- 
tional behavior for the streamwise pressure gradient 

/(,)— ~A 5 ln^ 7M -^j ) A 5 = (A9) 

As a conclusion, the two functional forms obtained assuming two-dimensionality of the flow 
i7j) and (|A9|) are not consistent with each other unless the function f(s) is a constant. 
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APPENDIX B: LEADING ORDER VELOCITY FIELD FOR ASYMMETRIC WALL 
PROFILE 



We present in this section the solution for the leading-order lubrication velocity field in the 
case of a planar channel of general shape; the calculation presented here contains that presented 
in section HirO as a special case. We use a cartesian coordinate system, as illustrated in Figured 
All notations are the same as in section lll Al with the difference that the channel is not right/left 
symmetric and we denote its boundaries by the equations y = hfi(x/X) and y = —hf2(x/\), 
with fi(x) > and f2{x) > 0. Using the same nondimensionalizations and regular expansion 
as in section III Al and defining the convenient notations fi = f\(x) and f2 — f2(x), we obtain 
that the leading-order axial velocity field is given by 

( \ ^Tul \ t 77 ( \ o dp (-ir / coshfc T1 (y + ^) \ 

Uo(.x,y,z) = > U n (x,y) cos k n z, U n (x,y) = 2- — — <^ , u + h , 1 > , (Bl) 

^ dx kl [ cosh/c n J 

and a pressure gradient obtained by enforcing the flow rate condition 

2/ [ fl u dydz = l, (B2) 
Jo J-h 

which leads to 

dp,, 3 ( 12 „ / tanh fr,*±&) \ V 1 

^^WTW)\TTT^{ « ) ) ' m 

Concerning the other two velocity components, the arguments presented in section IIIIBI also 
apply to the configuration considered here and the equations to solve for u and wo are equations 
f!38aj) and (|37|h along with the no-slip boundary conditions. Using the top/bottom symmetries 
in the velocity field 

v (x,y,-z) =v (x,y,z), w Q (x,y, -z) = -w Q (x,y, z) (B4) 
we look for for a solution of (|38a|) under the form of a double Fourier series in y and z, 

''o(-'-- !)■ '~ ) = ^2 An ( x > y) cos knZ + Bm ( x > z ) sin ( ^ 2y+ J 2 — —j (B5) 

n>0 m>0 ^ ' 



p>0 
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where k n = (n + l/2)7r, l m = run, k p = (p + l/2)ir and 



A n (x,y) = a n (x)P n (x,y) + b n (x)Q n (x,y), (B6a) 
B m (x,z) = c m (x)S rn (x,z), (B6b) 
C p (x,z) = d p (x)T p (x, z), (B6c) 



with 



P n (x, y) = (fx + f 2 ) sinh k n (f 1 - y) - [ft - y) cosh k n (fi - y) tanh k n (f 1 + f 2 ), (B7a) 

Q n (x,y) = (fi -y) sinh A; n (/i - y) - (fx - y) cosh k n {h - y) tanh k n (f 1 + f 2 ), (B7b) 

sinh ' 2/ : 



2£ m Z \ HUUi l/i+Ja 



S^x,^ = cosh ( - z ? v , (B7c) 

J1 + /2/ tanh I 



/1+/2 

2/-V- \ sinh fe 



A + /2 



T p (x,z) = cosh p - z V (B7d) 



tanh ' -2k 



,/l+/2 

Note that the solution in ()B5|) satisfies the no-slip boundary condition on all four walls. The 
integration of the continuity equation ()37|) leads to the general form of the vertical velocity 
component 



(B8) 



+ ^ £ m c m (x)Y m (x, z) cos 



m>0 



^m(2y + / 2 -/i) 
A + /a 



+ V" /cpc/p(x)Zp(x, 2) sin ( !?E2y+A — 111 

&o \ h + h 



with 



z cosh ' 2 ' : 



YJ.r. : ) - { ( 1 + t v I sinh ( -M^f J, (B9a) 

'jOtaiihfA-W V/1 + /2/ tanhte 



/1+/2 7 V/1+/2 7 / V-fr+'ft 

I I / , I \ . / 2/,'„: \ ^cosh(|k_ 



z ^> " = 5 1 + sua smh tei " > ' (B9b) 



The sets of unknown functions (a n ), (b n ), (c n ) and (d n ) are finally determined by enforcing the 
no-slip boundary condition for wq in (|B8j) on the four walls z = ±1, y — fi and y = —f 2 , and 
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we obtain a set of four infinite linear algebraic equations 

n>0 m>0 p>0 

b n (x) = £(*)+£52>(s)a i (s) + £^ (n>0), 

i>0 m>0 p>0 (BIO) 

c m (x) = c^(x) + y^(7^(x)ai(x) + J2C^l(x)b n (x)+J2C^(x)d p (x) (m>0), 

j>0 n>0 p>0 

= ^(x) + ^^(xKOr) + ^ D%>(x)b n (x) + ]T Dg(x)c m (x) (p > 0), 

i>0 n>0 m>0 



with 



/<9P \ _1 dv 

= ~ \di {xJi{x)) ) ^ (x ' /i(;r)) ' (Blla) 

W = - (^(*» 1 ( Bllb ) 
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and 



.4 



(2) 



(3) 



.4 



A (4) 

A ip 



(1) 



B 



^ 'mi 



(7(4) 



pi 



^pn 



£,(3) 

pm 



X 



X) 



X) 



X 



X) 



X) 



X 



X 



X 



X 



-5„ 



dPi 



dy 

2{-l) m k t i 



(x,h(x)) 
dp 



dy 



x,h{x)), 



2( l) p kik p 
'dQ n 



dy 
dy 



(xji(x)) 



-1 rl 



Y m (x, z) s\n(k n z) dz, 



-1 r l 



dy 



%,-h{x)) 



-i 



dPn 

dy 



Z p (x, z) sm(k n z) dz, 



[x,-f 2 (x)), 



x) — 2( 1) k n t Ti 



2( l) p ~*~ k n kp 
2(-l) 



dQr 

dy 
dQ n 



-1 rl 



dy 



X, ~h{x)) 
X, ~f2(x)) 

h dp 



Y m (x,z) sm(k n z) dz, 



-i r i 



Z p (x, z) sm.{k n z) dz, 



ki£ m Y m (x, l)(/i + f 2 ) J- h dy 
2(-l) n [ h dQ n 



k n £ m Y m (x, l)(/i + f 2 ) J- h dy 



x, y) cos 
x, y) cos 



t (2y + / 2 -/i) 

/1 + /2 
'J2y + f 2 -h) 



dy, 



fi + h 



dy, 



-2k p Zp(x, 1) 



^ m y m (x,i)(/i4 


■/a) J 


2(-l) 4 




kik p Zp(x, l)(/i 


+-/ 2 ) 


2(-l) n 





/l 



k n k p Zp(x, l)(/i + / 2 ) 
2£ m i / m (x.l) ^" 
fc p Z p (:r, l)(/i + / 2 ) 



cos 



4>(2y + / a -/0 



/1 + /2 



sm 



k p (2y + f 2 -f 1 ) 



fl dP 



h 



dQ, 
dy 



x, y) sin 
y) sin 



kp{2y + f 2 -h) 
/i + / 2 
k P (2y + f 2 -h) 



h + h 



fi + h 



(B12a) 
(B12b) 
(B12c) 
(B12d) 
(B12e) 
(B12f) 
(B12g) 
(B12h) 
«12i) 
(B12j) 
(B12k) 



cos 



u^y+h-h) 



sm 



^(2, + / m dj . B121) 

h + h J 



h + h 

Similar numerical technique as the ones illustrated in section IIIIEI can be used to solve (|B10|) 
and obtain the leading-order velocities in the general case given by (jB5|) and (IB8I) . 
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FIGURES 
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FIG. 1: Generic view of a micro-channel constrained to remain between two parallel planes. The design 
of channel has two degrees of freedom: (1) the shape of its centerline and (2) the relative width of the 
channel around this centerline. 
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FIG. 3: Curved micro-channel of constant cross section and slowly varying planar curvature. 
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FIG. 4: Illustration of the leading-order three-dimensional flow in the straight planar channel of varying 
dimensionless cross-section given by f(x) = 1 + 0.7 sin x. Top: view the channel. Bottom: plots of 
the leading-order dimensionless cross-sectional velocity field {vq,wq) and axial vorticity u)q at three 
locations along the channel: x = 5.25, 6.15 and 7.05; (a): in-plane velocity plots (the velocities are 
normalized by their maximum in-plane values); (b): iso- values of the y-component vq of the velocity, 
from equation (|43|): (c): iso- values of the z-component of the velocity Wo, from equation (|47j): (d): 
iso- values of the ^-component of the vorticity ujq, from equation 1)51 [I. 
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FIG. 5: Three-dimensional leading-order streamlines in the planar channel of varying dimensionless 
cross-section given by f(x) = 1 + 0.7 sin x. The channel is the same as the one illustrated in Figure 
|l]and only the streamlines in the quadrant (y > 0, z > 0) are reported; those in the other quadrants 
can be found using the flow symmetries (|42|) . The dimensionless time step used for computation is 
0.025 and 35 initially evenly spaced streamlines are considered. Top: three-dimensional view of the 
streamlines between the location of minimum width x = 3tt/2 (squares, filled) and the location of 
maximum width x = 5tt/2 (diamonds); the channel boundary and centerplane are also displayed. 
Bottom: (a) projection of the streamlines onto the (y, z) plane; (b) projection of the streamlines onto 
the (x, z) plane. 
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FIG. 6: Illustration of the leading-order three-dimensional flow strength in the planar channel of 
varying dimensionless cross-section f(x) = 1 + 0.7 sin x. The channel is the same as the one illustrated 
in Figure 0] and only the velocities in the quadrant (y > 0, z > 0) are considered. Figures (a) and (b): 
maximum cross-sectional values of the three components of the leading-order dimensionless velocity 
along the channel uq (circles), vq (squares, filled) and wq (triangles); (a): regular scale, (b): semi- log 
scale; note that when vq and wo were found to be zero, which happens at each location along the 
channel where f'(x) = under the lubrication approximation, they were replaced by 10 -5 for the 
semi- log figures. Figures (c) and (d): same as in (a) and (b) for the quadrant-averaged velocities 
< uq >, < vq > and < wq >; (c): regular scale, (d): semi-log scale. 
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